
clear all
************************************
*Numeric Export disabled
local TEX ""
local TEXON 1 //Change "0" to "1" to enable exporting numeric results to .txt files
if !`TEXON' local TEX "*"
*Also need to install `texresults2
*ssc install texresults2
************************************

***How much of difference in covariance is due to sampling error?

************************
*Set working directory
************************

program randomerrorcov
use "Enns_SociologicalScience_Reproduction\DS1_0.dta", clear
*generate series to mimic
*low, middle, and high income groups

quietly: sum pred90_sw
local high_mean = `r(mean)'
local high_min = `r(min)'
local high_max = `r(max)'

quietly: sum pred10_sw
local low_mean = `r(mean)'
local low_min = `r(min)'
local low_max = `r(max)'

quietly: sum pred50_sw
local mid_mean = `r(mean)'
local mid_min = `r(min)'
local mid_max = `r(max)'

//Create a matrix with the mean values from above.
matrix mean = (`high_mean', `mid_mean', `low_mean')

//Generate coveriance matrix
quietly: cor pred90_sw pred50_sw pred10_sw, cov
mat cov = r(C)

*Create 161 observations to represent 161 pairs in Gilens analysis of error covariance
*will drop values outside observed min and max, so need more than 161 here.
drawnorm high mid low, n(168) means(mean) cov(cov) clear

drop if high<`high_min' | high>`high_max' | low<`low_min' | low>`low_max'| mid<`mid_min' | mid>`mid_max'

*generate two versions based on above DGP
*the only difference is estimated sampling error
*.098 based on MOE if n=100, would be used to calculate deciles
*rnormal(0, 0.05) ///normally distributed w/ mean 0 and 95% of values +/- 0.098
//These results most likely an underestimate, because the income preferences are imputed, which would produce additional error.


gen high1 = high+rnormal(0, 0.05)
replace high1 = 0 if high1<0 //make sure no values exceed 0 or 1
replace high1 = 1 if high1>1 //make sure no values exceed 0 or 1
gen high2 = high+rnormal(0, 0.05)
replace high2 = 0 if high2<0
replace high2 = 1 if high2>1
gen mid1 = mid+rnormal(0, 0.05)
gen mid2 = mid+rnormal(0, 0.05)
replace mid1 = 0 if mid1<0
replace mid1 = 1 if mid1>1
replace mid2 = 0 if mid2<0
replace mid2 = 1 if mid2>1
gen low1 = low+rnormal(0, 0.05)
gen low2 = low+rnormal(0, 0.05)
replace low1 = 0 if low1<0
replace low1 = 1 if low1>1
replace low2 = 0 if low2<0
replace low2 = 1 if low2>1
end

*https://www.random.org/
set seed 60004172
tempfile randomcoverror
tempname temp1
postfile `temp1' simnum highmidcoverr highlowcoverr midlowcoverr n ///
			using "`randomcoverror'", replace

forval x = 1/1000 {
*run program that simulates series
randomerrorcov
	*generate covariance matrix for each version of question
cor high1 mid1 low1, cov
mat cov1 = r(C)
local highmid1 = cov1[2,1]
local highlow1 = cov1[3,1]
local midlow1 = cov1[3,2]
cor high2 mid2 low2, cov
mat cov2 = r(C)
local highmid2 = cov2[2,1]
local highlow2 = cov2[3,1]
local midlow2 = cov2[3,2]
*identify proportion difference

local coverrhighmid = `highmid1'/`highmid2'
local coverrhighlow = `highlow1'/`highlow2'
local coverrmidlow = `midlow1'/`midlow2'

local propcoverrhighmid = abs(1-`coverrhighmid')
local propcoverrhighlow = abs(1-`coverrhighlow')
local propcoverrmidlow = abs(1-`coverrmidlow')
count
local count = r(N)
post `temp1' (`x') (`propcoverrhighmid') (`propcoverrhighlow') (`propcoverrmidlow') (`count')
}


postclose `temp1' //write results to temp dataset 
use "`randomcoverror'", clear

centile highlowcoverr, centile(2.5 50 97.5)
local prophighlowcovlow = 100*`r(c_1)'/.19 
local prophighlowcovmed = 100*`r(c_2)'/.19 
local prophighlowcovup = 100*`r(c_3)'/.19 

centile highmidcoverr, centile(2.5 50 97.5)
local prophighmidcovlow = 100*`r(c_1)'/.17
local prophighmidcovmed = 100*`r(c_2)'/.17
local prophighmidcovup = 100*`r(c_3)'/.17

centile midlowcoverr, centile(2.5 50 97.5)
local propmidlowcovlow = 100*`r(c_1)'/.14
local propmidlowcovmed = 100*`r(c_2)'/.14
local propmidlowcovup = 100*`r(c_3)'/.14



`TEX'texresults2 using TxtFiles_NumericalResults\randommoe.txt, texmacro(prophighlowcovlow) result(`prophighlowcovlow') round(1) replace //
`TEX'texresults2 using TxtFiles_NumericalResults\randommoe.txt, texmacro(prophighlowcovmed) result(`prophighlowcovmed') round(1) append //
`TEX'texresults2 using TxtFiles_NumericalResults\randommoe.txt, texmacro(prophighlowcovup) result(`prophighlowcovup') round(1) append //
`TEX'texresults2 using TxtFiles_NumericalResults\randommoe.txt, texmacro(prophighmidcovlow) result(`prophighmidcovlow') round(1) append //
`TEX'texresults2 using TxtFiles_NumericalResults\randommoe.txt, texmacro(prophighmidcovmed) result(`prophighmidcovmed') round(1) append //
`TEX'texresults2 using TxtFiles_NumericalResults\randommoe.txt, texmacro(prophighmidcovup) result(`prophighmidcovup') round(1) append //
`TEX'texresults2 using TxtFiles_NumericalResults\randommoe.txt, texmacro(propmidlowcovlow) result(`propmidlowcovlow') round(1) append //
`TEX'texresults2 using TxtFiles_NumericalResults\randommoe.txt, texmacro(propmidlowcovmed) result(`propmidlowcovmed') round(1) append //
`TEX'texresults2 using TxtFiles_NumericalResults\randommoe.txt, texmacro(propmidlowcovup) result(`propmidlowcovup') round(1) append 

